{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## 第9章-EM算法及推广-导读\n",
    "\n",
    "&emsp;&emsp;从第9章开始，我们关注的点就和前面的章节不一样了，在第2-8章是分类问题，都属于监督学习，第9章EM算法是非监督学习，本章的推导比较多，在导读部分，重点介绍EM算法的应用、处理问题的特点、与之前算法的区别、EM算法的解决问题流程以及EM算法的性质和简单的变形。  \n",
    "\n",
    "### EM算法的引入\n",
    "&emsp;&emsp;EM算法是用来估计概率分布的，该概率分布的数据是有缺失的，假如数据没有缺失，$X \\sim N(\\mu, \\sigma^2)$，其中$\\mu,\\sigma^2$是未知的，用一组数据$(x_1,x_2,\\cdots,x_N)$去估计$\\mu,\\sigma^2$，估计方法一般为极大似然估计和贝叶斯估计，以上是正常估计密度函数所用的方法。  \n",
    "&emsp;&emsp;估计随机变量的密度函数属于无监督学习，但是在EM算法中，观测数据是有缺失的，书中有一个例子“三硬币模型”。  \n",
    "\n",
    "#### EM算法\n",
    "##### 三硬币模型\n",
    "&emsp;&emsp;随机变量$Z \\sim b(1,\\pi)$表示观测不到的数据，以概率$\\pi$取值为1，以概率$1-\\pi$取值为0，第1次实验得到$z_1$，根据该值得到观测值$y_1$，如果$z_1=1$，则$y_1 \\sim b(1,p)$，如果$z_1=0$，则$y_1 \\sim b(1, q)$，已知$y_1$，并不知道$z_1$取0或1，只知道$z_1 \\sim b(1,\\pi)$，$\\pi,p,q$未知。  \n",
    "&emsp;&emsp;根据上面的规则可得到$(z_2, y_2), \\cdots, (z_N, y_N)$，令$\\theta=(\\pi,p,q)$，将$(z,y)$称为完全数据，观测数据$y$称为不完全数据，所以能看到的就是不完全数据，对应完全数据有概率分布$P(y,z)=P(z)P(y|z)$，需要估计的是$P(y,z)$，也就是$\\pi,p,q$。一般使用的是极大似然估计，需要写出观测值：$$P(y|\\theta)=\\sum_z P(y,z|\\theta)=\\pi p^y (1-p)^{1-y} + (1-\\pi) q^y (1-q)^{1-y}$$\n",
    "所有的观测值$Y$为$$\\prod_{i=1}^N P(y_i|\\theta)=\\prod_{i=1}^N \\big[\\sum_z P(y_i,z|\\theta) \\big]=\\prod_{i=1}^N \\big[\\pi p^{y_i} (1-p)^{1-y_i} + (1-\\pi) q^{y_i} (1-q)^{1-y_i} \\big]$$\n",
    "&emsp;&emsp;极大似然估计就是极大化似然函数$\\displaystyle \\prod_{i=1}^N P(y_i|\\theta)$，为了简单一点可以取对数，但是后面的形式非常复杂，直接求导等于0非常困难，这就需要EM算法，用迭代的方式求解。\n",
    "\n",
    "##### EM算法的解释\n",
    "&emsp;&emsp;EM算法并不是直接处理似然函数的，为了简化这个问题，将最大化观测值的似然函数转换为最大化完全数据的似然函数，书中给出了两个角度的解释，第1种角度的解释在9.1节，第2种在9.4节。  \n",
    "完全数据的似然函数：$$\\ln \\prod_{i=1}^N P(y_i,z_i|\\theta)$$\n",
    "&emsp;&emsp;但是现在有一个问题：$z_i$未知。一旦给出$y,z$分布的形式，就可以将$y,z$的联合概率密度写成关于$\\theta$的函数，每一个$y_i$的值是已知的，但对于$z_i$的只是未知的，要如何处理$z_i$，就可以利用EM算法中的E步解决：将$z_i$有关的项，用期望代替$z_i \\rightarrow E(z)$，对应不同的$z$，期望$E(z)$都是相同的，求期望的时候需要知道密度函数（随机变量的分布），这个分布由估计量$\\theta,y_i$决定。给定$y_i,\\theta$求$z_i$的分布可以用贝叶斯公式，$y_i$已知，但是$\\theta$是需要估计的，所以这个时候采用迭代的方法。  \n",
    "**E步：** 首先初始化$\\theta^{(0)}$，然后进行迭代，第$i$步记为$\\theta^{(i)}$，这样就可以求出期望$E(z_i)$，将期望代入似然函数中。  \n",
    "**M步：** 最大化似然函数，求解$\\displaystyle \\theta^{(i+1)}=\\mathop{\\arg \\max} \\limits_{\\theta} \\ln \\prod_{i=1}^N P(y_i,E(z)|\\theta)$  \n",
    "重复E步，求期望$E(z)$，再继续M步，不断迭代。  \n",
    "**总结：** 在EM算法中，将极大化观测数据的似然函数转换为极大化完全数据的似然函数，在该过程中，隐变量$z_i$的值是未知的，采用给定参数的条件下，$E(z)$代替每一个$z_i$的值，得到用期望代替之后的似然函数，最大化该似然函数，一步一步对$\\theta$进行迭代。  \n",
    "\n",
    "##### EM算法\n",
    "输入：观测变量数据$Y$，隐变量数据$Z$，联合分布$P(Y,Z|\\theta)$，条件分布$P(Z|Y,\\theta)$  \n",
    "输出：模型参数$\\theta$  \n",
    "(1)选择参数的初值$\\theta^{(0)}$，开始迭代  \n",
    "(2)**E步：**记$\\theta^{(i)}$为第$i$次迭代参数$\\theta$的估计值，在第$i+1$次迭代的E步，计算$$\\begin{aligned} Q(\\theta,\\theta^{(i)})\n",
    "=& E_Z\\big[\\ln P(Y,Z|\\theta) | Y, \\theta^{(i)}\\big] \\\\\n",
    "=& \\sum_Z \\ln P(Y,Z|\\theta)P(Z|Y,\\theta^{(i)})\n",
    "\\end{aligned}$$这里，$P(Z|Y,\\theta^{(i)})$是在给定观测数据$Y$和当前的参数估计$\\theta^{(i)}$下隐变量数据$Z$的条件概率分布  \n",
    "(3)**M步：**求使得$Q(\\theta,\\theta^{(i)})$极大化的$\\theta$，确定第$i+1$次迭代的参数估计值$\\theta^{(i+1)}$ $$\\theta^{(i+1)}=\\mathop{\\arg  \\max} \\limits_{\\theta} Q(\\theta,\\theta^{(i)})$$ \n",
    "(4)重复第(2)步和第(3)步，直到收敛（收敛条件：$\\theta^{(i)}$和$\\theta^{(i+1)}$很接近，或者是$Q(\\theta^{(i+1)},\\theta^{(i)})$和$Q(\\theta^{(i)},\\theta^{(i-1)})$很接近）。  \n",
    "函数$Q(\\theta,\\theta^{(i)})$是EM算法的核心，称为$Q$函数。  \n",
    "\n",
    "#### EM算法的导出\n",
    "<br/><center><img style=\"border-radius: 0.3125em;box-shadow: 0 2px 4px 0 rgba(34,36,38,.12),0 2px 10px 0 rgba(34,36,38,.08);\" src=\"../../../PhaseFour/Note/image/9-1-Explanation-of-EM-algorithm.png\"><br><div style=\"color:orange; border-bottom: 1px solid #d9d9d9;display: inline-block;color: #000;padding: 2px;\">图9-1 EM算法的解释</div></center>  \n",
    "\n",
    "&emsp;&emsp;图9-1中实线为观测数据的似然函数$L(\\theta)=P(y|\\theta)$，极大似然函数寻找该曲线的最大值，将这个点对应的$\\hat{\\theta}$作为$\\theta$的估计值，$Q(\\theta,\\theta^{(i)})$函数：当给定$\\theta^{(i)}$时，该函数是一个关于$\\theta$的函数，图中是用$B(\\theta,\\theta^{(i)})$来代替的，$B(\\theta,\\theta^{(i)})$和$Q(\\theta,\\theta^{(i)})$取极值的点是一样的，$B(\\theta,\\theta^{(i)})$是$L(\\theta)$的下界，求解$\\displaystyle \\hat{\\theta}^{(i)} = \\mathop{\\arg \\max} \\limits_{\\theta} B(\\theta,\\theta^{(i)})$，这样就可以求解下一个点$\\theta^{(i+1)}$使得函数$B$最大化，求解${\\hat{\\theta}}^{(i+1)}$，不断迭代。  \n",
    "&emsp;&emsp;从图中看，$L(\\theta)$并不是一个凸函数，如果求极大化，该函数并不是一个凹函数，也就是说，用这种迭代的方法求$L(\\theta)$极大值时，当选取的初值不同，可能会收敛到不同的极大值处，所以EM算法得到的结果和初值的选择有关，初值选择不同，可能会得到不同的极大值点。  \n",
    "\n",
    "### EM算法的收敛性\n",
    "\n",
    "> **定理9.1** 设$P(Y|\\theta)$为观测数据的似然函数，$\\theta^{(i)}(i=1,2,\\cdots)$为EM算法得到的参数估计序列，$P(Y|\\theta^{(i)})(i=1,2,\\cdots)$为对应的似然函数序列，则$P(Y|\\theta^{(i)})$是单调递增的，即$$P(Y|\\theta^{(i+1)}) \\geqslant P(Y|\\theta^{(i)})$$**定理9.2** 设$L(\\theta)=\\log P(Y|\\theta)$为观测数据的对数似然函数，$\\theta^{(i)}(i=1,2,\\cdots)$为EM算法得到的参数估计序列，$L(\\theta^{(i)})(i=1,2,\\cdots)$为对应的对数似然函数序列。  \n",
    "(1)如果$P(Y|\\theta)$有上界，则$L(\\theta^{(i)})=\\log P(Y|\\theta^{(i)})$收敛到某一值$L^*$  \n",
    "(2)在函数$Q(\\theta,\\theta')$与$L(\\theta)$满足一定条件下，由EM算法得到的参数估计序列$\\theta^{(i)}$的收敛值$\\theta^*$是$L(\\theta)$的稳定点。\n",
    "\n",
    "&emsp;&emsp;本节给出了两个定理：第一个定理，用EM算法求解得到的$\\theta$，每次更新$\\theta^{(i)}$时，都会使得观测数据的似然函数变大，这个性质可以保证求解最大化的$\\theta$值。将最大化观测数据的似然函数转换为最大化完全数据的似然函数，但是通过最大化完全数据的似然函数得到的$\\theta$的更新值，会使得观测数据的似然函数越来越大。第二个定理，对EM算法收敛性的一个说明，EM算法会收敛到一个局部的最大值。\n",
    "\n",
    "### EM算法在高斯混合模型学习中的应用\n",
    "&emsp;&emsp;高斯混合模型是一个非常常见的模型，其密度估计问题也是一个比较复杂的问题，但是可以用EM算法解决。  \n",
    "\n",
    "#### 高斯混合模型\n",
    "定义9.2（高斯混合模型） 高斯混合模型是指具有如下形式的概率分布模型：$$P(y|\\theta)=\\sum_{k=1}^K \\alpha_k \\phi(y|\\theta_k)$$其中，$\\alpha_k$是系数，$\\displaystyle \\alpha_k \\geqslant 0, \\sum_{k=1}^K \\alpha_k = 1$，$\\phi(y|\\theta_k)$是高斯分布密度，$\\theta_k=(\\mu_k, \\sigma_k^2)$，$$\\phi(y|\\theta)=\\frac{1}{\\sqrt{2 \\pi} \\sigma_k} \\exp \\left( -\\frac{(y-\\mu_k)^2}{2\\sigma_k^2} \\right)$$称为第$k$个分模型。  \n",
    "&emsp;&emsp;首先介绍高斯混合模型，只考虑简单的一维随机变量$y$，高斯分布就是正态分布，$y \\sim N(\\mu, \\sigma^2)$，给定$y$的观测值，就可以很容易求出$\\mu$和$\\sigma^2$，但是目前$y$不是来自高斯分布，而是有一定概率的来自于两个不同的高斯分布$N(\\mu_1, \\sigma_1^2)$和$N(\\mu_2, \\sigma_2^2)$，这个就是两个高斯分布的混合，并不知道$y$来自于哪一个高斯分布，这里涉及到了**隐变量**，这里所举的例子是两个高斯分布的混合，对于一般的形式如定义9.2给出。  \n",
    "&emsp;&emsp;$\\theta$包含了两个部分，第1部分是取值属于哪个高斯分布，模型中一共有$K$个高斯分布，另一部分是多项分布，一共有$K$个取值，只可能取到其中的一个值，对于不同的取值，都对应了一个概率，这个概率的和等于1，取不同的$k$个概率，对应的就是$\\alpha_k$，也就是以$\\alpha_1$的概率取了第1个高斯分布，以$\\alpha_2$的概率取了第2个高斯分布，以此类推。$\\phi(y|\\theta_k)$表示第$k$个高斯分布，其参数对应的$\\mu,\\sigma^2$为$\\theta_k$，在给定参数$\\theta_k$时，$P(y|\\theta)$为所求。  \n",
    "&emsp;&emsp;三个硬币的例子中，其模型为$\\displaystyle P(y|\\theta)=\\sum_Z p(y,z|\\theta)=\\sum_Z p(z|\\theta) P(y|z,\\theta)$，对应于高斯混合模型中，$p(z|\\theta)$是一个多项分布，一共有$K$项，且每一项的概率都为$\\alpha_k$，$P(y|z,\\theta)$就是确定了某一个$k$的取值之后的高斯分布，把之前的这两个公式代入之后求和，推导之后得到的公式为$\\displaystyle P(y|\\theta)=\\sum_{k=1}^K \\alpha_k \\phi(y|\\theta_k)$  \n",
    "\n",
    "#### 高斯混合模型的参数估计\n",
    "&emsp;&emsp;对于高斯混合模型的参数估计，书中给出了详细的推导过程，这里就简单的介绍一下，在EM算法中，极大化完全数据的似然函数$P(y,z|\\theta)$，在完全数据中，$z$是观测不到的，所以在每个E步中都给定当前的$\\theta$值和$y$，得到期望$E(z|y,\\theta^{(i)})$，用期望代替$z$，最后最大化转换之后的似然函数，得到$\\theta$的下一个估计值，这就是EM算法在高斯混合模型中的一个流程。  \n",
    "&emsp;&emsp;在该模型中，$z$是一个分类的变量，取值为$1, 2, \\cdots, K$，书中的推导过程是用了一个向量$\\gamma$表示$z$，如果$z=1$，则$\\gamma=(1,0,0,\\cdots,0)$，如果$z=2$，则$\\gamma=(0,1,0,\\cdots,0)$，这个相当于One-hot，也就是说$z$是第$i$个高斯分布，在$\\gamma$的第$i$个分量为1，其他分量都为0，以上是书中对隐变量的处理。  \n",
    "\n",
    "### EM算法的推广\n",
    "&emsp;&emsp;本节从另一个角度解释了EM算法，EM算法可以解释为$F$函数的极大-极大算法。  \n",
    "\n",
    "> **算法9.4（GEM算法2）**  \n",
    "输入：观测数据，$Q$函数  \n",
    "输出：模型参数  \n",
    "(1)初始化参数$\\theta^{(0)}$，开始迭代  \n",
    "(2)第$i+1$次迭代，第1步：记$\\theta^{(i)}$为参数$\\theta$的估计值，计算$$\\begin{aligned} Q(\\theta,\\theta^{(i)})\n",
    "=& E_Z\\big[ \\log P(Y,Z|\\theta)|Y,\\theta^{(i)} \\big] \\\\\n",
    "=& \\sum_Z P(Z|Y,\\theta^{(i)}) \\log P(Y,Z|\\theta)\n",
    "\\end{aligned}$$\n",
    "(3)第2步：求$\\theta^{(i+1)}$使得$$Q(\\theta^{(i+1)},\\theta^{(i)}) > Q(\\theta^{(i)},\\theta^{(i)})$$\n",
    "(4)重复(2)和(3)，直到收敛  \n",
    "\n",
    "> **算法9.5（GEM算法3）**  \n",
    "输入：观测数据，$Q$函数  \n",
    "输出：模型参数  \n",
    "(1)初始化参数$\\theta^{(0)}=(\\theta^{(0)}_1,\\theta^{(0)}_2,\\cdots,\\theta^{(0)}_d)$，开始迭代  \n",
    "(2)第$i+1$次迭代，第1步：记$\\theta^{(i)}=(\\theta^{(i)}_1,\\theta^{(i)}_2,\\cdots,\\theta^{(i)}_d)$为参数$\\theta=(\\theta_1,\\theta_2,\\cdots,\\theta_d)$的估计值，计算$$\\begin{aligned} Q(\\theta,\\theta^{(i)})\n",
    "=& E_Z\\big[ \\log P(Y,Z|\\theta)|Y,\\theta^{(i)} \\big] \\\\\n",
    "=& \\sum_Z P(Z|Y,\\theta^{(i)}) \\log P(Y,Z|\\theta)\n",
    "\\end{aligned}$$\n",
    "(3)第2步：进行$d$次条件极大化：  \n",
    "首先，在$\\theta^{(i)}_2,\\theta^{(i)}_3,\\cdots,\\theta^{(i)}_k$保持不变的条件下求使得$Q(\\theta,\\theta^{(i)})$达到极大的$\\theta^{(i+1)}_1$；  \n",
    "然后，在$\\theta_1=\\theta^{(i+1)}_1,\\theta_j=\\theta^{(j)}_j,j=3,4,\\cdots,k$的条件下求使$Q(\\theta,\\theta^{(i)})$达到极大的$\\theta^{(i+1)}$；  \n",
    "如此继续，经过$d$次条件极大化，得到$\\theta^{(i+1)}=(\\theta^{(i+1)}_1,\\theta^{(i+1)}_2,\\cdots, \\theta^{(i+1)}_d)$使得$$Q(\\theta^{(i+1)},\\theta^{(i)}) > Q(\\theta^{(i)},\\theta^{(i)})$$\n",
    "(4)重复(2)和(3)，直到收敛\n",
    "\n",
    "&emsp;&emsp;算法9.3和EM算法是完全一致的，只不过将$Q$函数写成了$F$函数形式。算法9.4中第(2)步和E步是一致的，第(3)步和M步是一致的，在EM算法中的M步是求解使得$Q$函数最大的$\\theta$，在算法9.4中，要求就放宽了一些，求解使得$Q$函数有提升的$\\theta$，这就是EM算法的第1个变形。对于EM算法的第2个变形（算法9.5），E步是不变的，在M步上做了改进，如果$\\theta$有很多不同的分量，分别求解使得$Q$函数最大的每一个分量，当求解第1个分量$\\theta_1$时，固定其他分量，以便求解分量$\\theta_1$，求解其他分量的方法相同，这样就将求解$\\theta$向量的最优解变成求解每一个分量的最优解，很多时候，这样的求解方法能大大简化计算。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.2"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
